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The dynamics and energetics of intergranular crack growth along a flat grain boundary in 
aluminum is studied by a molecular-dynamics simulation model for crack propagation under 
steady-state conditions. Using the ability of the molecular-dynamics simulation to identify atoms 
involved in different atomistic mechanisms, it was possible to identify the energy contribution of 
different processes taking place during crack growth. The energy contributions were divided as: 
elastic energy - defined as the potential energy of the atoms in fee crystallographic state; and 
plastically stored energy - the energy of stacking faults and twin boundaries; grain-boundary and 
surface energy. In addition, monitoring the amount of heat exchange with the molecular- 
dynamics thermostat gives the energy dissipated as heat in the system. The energetic analysis 
indicates that the majority of energy in a fast growing crack is dissipated as heat. This dissipation 
increases linearly at low speed, and faster than linear at speeds approaching 1/3 the Rayleigh 
wave speed when the crack tip becomes dynamically unstable producing periodic dislocation 
bursts until the crack is blunted. 
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1. Introduction 

Intergranular fracture is the governing process for mechanical failure in nanocrystalline metals 
[1]. Intergranular crack propagation is strongly influenced by the specifics of the grain boundary 
(GB) interface, the structure of which still poses many unresolved issues [2], The GB volume at 
the front of the approaching crack tip undergoes structural changes, which, in turn, affect the 
crack propagation. The presence of crystals with two different crystallographic orientations, 
containing slip systems oriented differently to the crack plane on both sides of the crack, is also 
another complicating factor which has to be taken into account. The presence of the GB also 
affects the energetics and dynamics of the intergranular crack propagation. While there are a 
large number of studies on the structural processes at the crack tip (some of them, related to 
intergranular cracks will be discussed below), little is known on the energy transformations 
accompanying these processes, which defines the crack growth resistance. 

The propagation of a crack involves a substantial amount of energy to be consumed or 
dissipated. The energy balance has a governing role in the process of fracture. In its simplest 
form, for a perfectly brittle material, the energy balance is expressed by Griffith’s equation [3] as 

G- 2y s (1) 

where G is the energy release rate associated with crack extension, and y s is the free energy of 
the newly created surface on each side of the crack. For ductile materials, Orowan [4] and Irwin 
[5] extended Eq. (1) to include the plastic energy per unit area y P i dissipated through various 
irreversible processes in the plastic zone surrounding the crack tip [6] 

G- 2y,+Y„. (2) 

In the case of intergranular fracture, the energy of the GB ygb has an additional contribution 
to the balance 
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G = 2y. s +y /;/ -y cb , (3) 

indicating that it is easier to grow a crack along a GB interface than through a perfect crystal. A 
systematic analysis on the relation between GB crystallography, GB energy, and GB cleavage 
properties was performed by Wolf [7] using molecular-dynamics (MD) simulations of copper 
and gold. In this study [7], a large representative set of symmetric tilt and twist GBs were 
considered according to two characteristic parameters, the twist and tilt angle (cp, 0), and their 
cleavage energies, E c /, defined by analogy with Eq. (3) as 

£ d= 2 Y,-Y GB (^ e ) ( 4 ) 

were estimated. As a result, a structure- energy correlation for all symmetric GBs in an fee metal 
was derived and related to the inter- granular fracture properties. The role of plasticity was not 
considered in ref. [7], but investigations on dislocation nucleation from the crack tip as a 
function of the GB crystallography and structure were performed in a series of subsequent works 
[8-10]. Intergranular crack induced plasticity in nanocrystalline metals was studied in [10-12]. 
The decohesion strength of a GB and its dependence on impurities was studied in [13-16]. The 
dependence of the fracture stress and strain of the GB interface on the atomic disorder due to the 
presence of vacancies and interstitials, has also been recently investigated [17]. 

One approach to analyze the plasticity contribution to crack resistance has been suggested by 
Stamp fl and Kolednik [18], The total plastically deformed volume is divided into several 
characteristic regions, according to the distance from the crack surface, and the energy 
contribution of each region to Ypi is determined. The idea is based on experimental findings that 
crack growth resistance often depends on the geometry and the size of the specimen, thus 
implying that different volume regions surrounding the crack contribute differently to G [18]. 
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A more fundamental physics based approach to study crack resistance is presented in this 
paper. The energy dissipation at the crack tip is determined through identifying the underlying 
physical process and quantifying individual contributions during crack propagation. A Molecular 
Dynamics (MD) simulation of a crack growing along a high-angle ^99 <11 0> tilt GB in 
aluminum is used to interrogate the atomic processes associated with the crack tip. MD 
simulations have the advantage that they can capture the dynamics of crack propagation at an 
atomic scale. By integrating the Newtonian equations of motion to calculate atomic trajectories, 
MD simulations can give invaluable information not only on the structural changes in the 
system, but also on the energetics and dynamics of the underlying processes. By using special 
techniques in identifying atoms in different structural states (specified further in the text), the 
energy spent for crack propagation can be separated into different categories according to its 
participation in certain physical processes. As a result, a complete picture of the energy balance 
of the crack propagation, revealing a variety of energy transformations can be obtained. Such an 
analysis could help in better understanding the energetics of crack propagation in nanocrystalline 
fee metals and in other materials, where intergranular fracture is a dominant failure mechanism. 

This paper is constructed as follows: The simulation approach is described in Section 2. The 
crack growth dynamics is presented in Section 3. A detailed analysis of the energy dissipated in 
the process of crack propagation is given in Section 4. The main conclusions of this study are 
outlined in Section 5. 


2. The simulation approach 

The simulation approach used in this study is based on a MD simulation model of crack 
propagation under time-independent, or steady-state conditions through a flat grain boundary in 
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aluminum modeled by the interatomic potential of Mishin et al. [19]. The simulation is 
performed at a temperature of 100 K to suppress GB and surface diffusion and facilitate brittle 
fracture. The temperature is maintained by using a molecular-dynamics Nose-Hoover thermostat 
[20]. The MD simulation set up follows the known analytical model for steady-state crack 
propagation [21] under constant strain conditions and has the configuration schematically shown 
in Fig. 1. With periodic boundary conditions in all directions, the model represents an aluminum 
multilayer system of alternating sets of thick and thin crystalline layers separated by four flat 
GBs. The two broad layers, marked as “Crystal I” and “Crystal II”, form a bicrystalline system 
with a flat GB in the middle, through which the crack propagates. The crystallographic 
orientations of Crystal I and Crystal II are presented in Fig. 2 representing a snapshot of an 
enlarged part of the simulation system around the GB after relaxation. In the imposed coordinate 

system of the model, the orientation of Crystal I is: (x:[7 7 10]; y:[5 5 7]; z: [1 1 0]), 

and the orientation of Crystal II is: (x:[ 7 7 10];y:[5 5 7];z:[l 1 0]) (see Fig. 2). In this 
way, Crystal II is a mirror image of Crystal I relative to the crystallographic plane {5 5 7}, 

which becomes the plane of the GB between them. The GB thus formed is a <1 1 0> ]T99 
symmetric tilt GB. Its atomic structure in Al agrees well with that determined through 
transmission electron microscopy (TEM) investigations [22]. This is a high- angle grain boundary 
(tilt angle of 89.42°) with a large excess (i.e., above the perfect crystal) energy, y G B = 0.60±0.05 
J/m 2 , estimated here for a relaxed structure at T = 100 K. The high excess GB energy facilitates 
its decohesion [7]. The surface energy of the GB plane y s at 100 K is estimated in this study as y s 
= 0.952±0.010 J/m 2 and is in good agreement with experimental data for the Al surface energy. 

The two smaller layers, Absorbing Layer I and Absorbing Layer II, on both sides of the 
bicrystalline system (Fig. 1) have the same crystallographic orientations as Crystal II and Crystal 
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I, respectively. Consequently, the GBs formed by these layers are of the same crystallographic 
type as the GB between Crystal I and Crystal II. The purpose of these layers in the simulation is 
to serve as absorbers for the phonon waves [15] generated from the crack tips preventing their re- 
entrance back into the system through the imposed periodic boundary conditions, which may 
influence the crack propagation. The wave absorption is performed by applying a damping 
friction coefficient to the atoms in the absorbing layers. This damping effectively takes energy 
from the system and decreases the temperature of the layers. The applied Nose-Hoover 
thermostat [20] compensates for the absorbed energy by supplying an additional amount of 
energy to keep the average temperature constant at 100K. The sum of the energy taken by the 
absorbing layers and by the thermostat would make the total heat A Q released by the system in 
the process of crack growth. In addition to absorbing phonon waves, the GBs created between 
these layers and Crystals I and II act as sinks for dislocations that may be emitted and propagate 
from the crack tips during crack growth. In this way, the negative effect of the periodic boundary 
conditions in the v-direction creating periodic images of all crack-tip disturbances and 
influencing the crack propagation is suppressed. 

The interaction between individual atoms in the system is presented by a many-body 
embedded-atom method (EAM) potential of Mishin et al. [19] fitted to give the correct zero- 
temperature lattice constant, a 0 = 4.05 A, elastic constants, cohesive energy, vacancy formation 
energy, etc. Of particular importance for the simulation of fracture and dislocation plasticity is 
the close fit of the potential to the experimentally measured surface and stacking-fault energies. 
Potential-dependent parameters that will be needed in this study are the relaxed stable stacking- 
fault energy, y s f = 0.146 J/m 2 , the unstable stacking-fault energy, y us = 0.168 J/m 2 , as defined in 
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ref. [19], and the unstable twinning energy, y ut = 0.210±0.010 J/m 2 , estimated here according to 
the method described in ref. [23]. 

The system thickness h in the z-direction equals only 10(2 2 0) crystallographic planes 
(accounting for the symmetry of the fee lattice), or h = 10V2/2a o - 2.9 mn. This thickness is more 
than four times larger than the range of the interatomic potential, r c = 1.55a 0 = 0.628 mn [19], 
which is sufficient to prevent interference of the atoms with their periodic images and to preserve 
the local three-dimensional (3D) physics in the system. The small thickness in the z-direction 
allows the system size in the x- and y- directions to extend up to 21 <7 7 10>a o -120 mn, and 
24<5 5 7>a 0 - 97 mn, respectively (see Fig 1 and Fig. 2), while limiting the number of simulated 
atoms to 1,994,000, allowing the simulation to be carried out on a modest Beowulf cluster. 

As the dislocation activity is expected to be very important in this study, the choice of the [1 
1 0] direction as a columnar axis for this quasi-two dimensional set-up ensures that two slip 

planes (the (1 1 1) and (1 1 1) planes) with a total of six slip systems are presented in each 
crystal [24,25], The existence of two slip planes allows for a variety of dislocation interactions 
and cross slip events to take place [24,25] resembling closely a full 3D environment, the main 
constraint being that the dislocation lines of all possible dislocations have to be straight lines 
parallel to the columnar direction. The implication of this constraint affects mostly the stress 
barrier for dislocation nucleation making it slightly higher for a straight dislocation, compared to 
a dislocation loop in a fully 3D system [26]. While neglecting the effects of curvature is not 
expected to qualitatively change the fracture mechanism, it may affect the process quantitatively 
in terms of slightly decreasing peak stress* and work of decohesion. 

* 

The stress to nucleate a dislocation from a grain boundary in a nanocrystalline columnar simulated microstructure 
was found to be 2.3 GPa [24], while in a fully 3D microstructure of the same material it was found to be less than 
2.0 GPa [26]. 
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After equilibration at zero constant pressure using the Parrinello-Rahman constant-stress 
simulation [27], the system is loaded hydrostatically in tension, i.e., 


°xx = °yy=°zz =0 


(5) 


and is dynamically equilibrated at this constant stress. After establishing equilibrium between the 


strain in the system and the applied external stress, the system size in all three dimensions is 


fixed under constant strain conditions. 


The transition from a constant stress to a constant strain simulation transforms the volume 


fluctuations that are always present in a finite system under thermodynamic equilibrium into 
stress fluctuations. Thus, further equilibration at constant strain is necessary to smooth out these 
fluctuations. The simulation then proceeds under this constant strain - constant temperature 
condition. Although the simulations are carried out under constant strain, for convenience in 
presentation, the various analyses will reference the value of prestress that corresponds to a 
particular value of prestrain. 

The crack in the system is nucleated by screening (preventing) the atomic interactions 
between atoms at both sides of the GB plane between Crystal I and Crystal II along a region of 
length / 0 . The crack starts growing if / 0 is larger than the critical Griffith length L„ defined when 
the energy spent to create the upper and lower crack surface 2y s minus the energy gained by 
destroying the GB y G B is equal to the released strain energy -dU/dl, per length /, 

2y s - Ygb = -dU/dl. (6) 

An estimate of L g is made by calculating dU/dl as a function of / and a. This calculation is 
performed by using an anisotropic elastic finite element model of the elastic equivalent for the 
MD system, as described in ref. [28]. L g sets the minimum size to grow a crack in the system. In 
order to simulate a sufficient crack growth length, L g must be at least one order of magnitude 
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smaller than the system size. For this reason and to reflect the existing periodicity of the 
equilibrated GB structure induced by the lattice of the two joined crystals [22], the screened 
region length is set equal to the crystallographic period in the x-direction, or 
l 0 = (7 7 10}a 0 = 5.7 mn. The performed calculations for L g satisfying Eq. (6) give L g < l 0 for 
a > 3.5 GPa [28]. Thus, the values of the prestress in this study were set to a = 3.5, 3.75, 4.0, and 
4.25 GPa. Crack growth has been observed even at a = 3.5 GPa because the interaction range of 
the many body potential increases the effective length of the screened region as all the atoms 
within the interaction range of the screened atoms are also affected by the screening and make 
the initial crack length effectively larger than prescribed. 

These very high prestress values, imposed by the severe system size limitations in the MD 
simulation, are larger than the yield strength of the material and would have caused strong 
plasticity effects not related to the crack, such as spontaneous dislocation nucleation from the 
GBs [24-26]. Applying triaxial hydrostatic stress (Eq. (5)) eliminates these undesirable plasticity 
effects. 

The identification of various structural defects including dislocations, twins, stacking faults, 
etc., appearing around the growing crack is important for understanding the various deformation 
mechanisms in the system. A procedure for atom identification based on the atom’s coordination 
number and on the common-neighbor-analysis (CNA) technique [29, 30] is used. The technique 
makes it possible to identify atoms in fee and hep states. Layers of hep atoms in an fee lattice are 
formed at stacking faults and twin boundaries [31] and can be used successfully for visualizing 
the ongoing dislocation processes in fee crystals [32]. Atoms that are not identified in an fee or 
hep state are considered to be in a non- crystalline state and indicate the presence of GBs or 
dislocation cores. In addition, atoms that have lost more than 1/3 of their neighbors inside the 
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interaction range of the potential are considered surface atoms. Under this convention, the 
thermalized structure of a (110)299 symmetric tilt GB, shown in Fig. 2, appears as quasi- 
periodic, with a regular pattern of hep atoms immersed in a disordered layer with few distributed 
vacancies, identified by the appearance of few isolated surface atoms inside the GB region. As 
will be shown later in this paper, classifying atoms in this way presents a unique possibility to 
distinguish and quantify the various atomic processes occurring at the crack tip, and to identify 
their contribution to the work of decohesion for the GB interface. 

3. Dynamic crack growth along the 299 grain boundary 

The process of intergranular crack growth at 3.75 GPa prestress is monitored through the 
series of snapshots in Fig. 3. The identified hep atoms form the stacking faults and twin 
boundaries visualized as dark straight lines (© in Figs. 3(a-d)) passing through the interiors of 
the two crystals. The disordered atoms form the GB seen as a thick black line dividing the two 
crystals (® in Fig. 3(a)), the dislocation cores of the twinning dislocations, seen along the twin 
boundaries (® in Figs. 3(b-d)), and the cores of the slip dislocations found in the crystals’ 
interior (© in Fig. 3(i)). The wavy patterns in the crystal regions are Moire patterns produced by 
the overlapping of the elastically deformed atomic lattice with the regular pixel grid of the digital 
image. In a way, these patterns make visible the strain patterns and phonon waves in the crystal. 

The MD simulation is used to determine the mechanisms of formation and crack growth from 
initiation at t = 0 ps through t = 104 ps. As seen in Fig. 3, the mechanisms of crack growth in the 
-x direction and +x direction are entirely different. 

As early as 4 ps after initiation (t = 4 ps), a nanotwin starts to form at the crack tip 
propagating in the -x direction (© in Fig. 3(a)). The nanotwin grows (© in Figs. 3(b, c)) by 
continuous emission of a series of twinning dislocations from the crack tip (® in Figs. 3(b, c)) 


10 



Invited paper for the Special Issue of Journal of Materials Science 
Special Issue Title: “Nanostructured Materials and Mechanical Behavior” 

[31]. At t= 12 ps after the crack initiation, a second nano twin at the same crack tip, but on the 
opposite side of the GB is nucleated (see both identifiers © in Fig. 3(c)) and proceeds to grow 
symmetrically with the path of the first nanotwin (© and © in Figs. 3(c, d)). Both nanotwins 
continue to grow in a symmetric fashion as shown in the remaining snapshots (Figs. 3(e-i)). 

Conversely, the crack propagation in the +x direction takes place through a continuous 
process of void formation and void coalescence (© in Figs. 3(d-i)). Compared to the crack 
growth in the -x direction, where the crack tip is effectively blunted by the nanotwins, the crack 
growth in the +x direction is very rapid. In spite of its intensity, this growth has not produced any 
dislocation emission for 76 ps after the crack initiation (Figs. 3(a-g)). In the time interval 
between 76 and 104 ps (see Figs. 3(g) and 3(i)), the crack tip suddenly emits a burst of 
dislocations (© in Fig. 3(i)) and the growth in the +x direction stops. Further simulation shows 
no more crack propagation or growth of the nanotwins, indicating that the crack has reached 
equilibrium at a length close to 40 mn. This asymmetric crack growth in the +x and -x directions 
is due to the specific orientation of the {111} slip planes on both sides of the GB interface 
forming different inclination angles with the propagation directions of the two crack tips. The 
observation is in agreement with the Rice [33] and Tadmor-Hai [23] criteria for twinning vs. 
cleavage at the crack tip, as demonstrated in ref [28], 

The equilibrium states for the other three cases of prestress, a = 3.5, 4.0, and 4.25 GPa are 
shown in Figs. 4(a-c). At the lowest prestress, a = 3.5 GPa (Fig. 4(a)), the crack exhibited very 
little growth (from I 0 = 5.7 mn nucleation length to / ~ 10 mn equilibrium length). The 
equilibrium was reached due to the formation of the two nanotwins at the tip (©in Fig. 4(a)), 
propagating in the -x direction, which plastically relieved enough strain energy to arrest the 
crack growth. At a = 4.0 GPa (Fig. 4(b)), the growth process is very similar to the case of a = 
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3.75 GPa, discussed earlier. In the -x direction, the crack is blunted by defonnation twins, 
growing symmetrically in the two crystals, while in the +x direction, the crack is arrested by a 
dislocation burst. Increasing the prestress to a = 4.25 GPa resulted in a substantial increase of the 
crack length (Fig. 4(c)), accompanied by more intensive plastic processes such as enhanced 
deformation twinning (© in Fig. 4(c)) and enhanced dislocation activity (CD and © in Fig. 4(c)). 
A more detailed description and analysis of the atomistic processes taking place at the two tips of 
the growing crack is given in ref. [28]. The present discussion will focus on the dynamics of the 
crack growth and the energy transformations at the crack tips which define the crack growth 
resistance. 

The speed of crack growth can be estimated in the MD simulation by determining the rate of 
increase of the crack surface area. This can be done by using the identification procedure for 
“surface” atoms as described in Section 2. On average, 1 nm 2 of a flat {5 5 7} surface contains 
16 surface atoms. The crack free surface S is estimated by counting the number of surface atoms, 
and when divided by 2 h (see Fig. 1, the prefactor 2 accounts for the two crack surfaces), 
conveniently gives an effective crack length /. 

The dependence of / with time, given in Fig. 5, shows distinctively different behavior for the 
four prestress values used in the study. At the lowest prestress, a = 3.5 GPa, the crack shows a 
slow, but steady growth at the speed of 25 m/s for almost 150 ps and then arrests (Fig. 5). The 
reason for this arrest is the formation of two symmetrical nanotwins at the crack tip propagating 
in the —x direction (Fig. 4(a)), in the same way as in the a = 3.75 GPa case, discussed previously 
(Fig. 3). These nanotwins relieve the stress around the crack [28] and prevent its further growth. 
The crack reaches an elasto -plastic equilibrium with the surrounding material at a length of 10 
mn (Fig. 4(a)). 
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At a slightly higher prestress of a = 3.75 GPa, the average speed of crack growth increases 
dramatically from 25 to 511 m/s (see the corresponding curve in Fig. 5). This fast growth takes 
place in the +x direction only, while the propagation in the —x direction is blunted by the active 
formation of the two symmetric nanotwins (© in Fig. 3). After reaching a length of 35 mn, the 
crack growth in the +x direction is arrested, as shown in Fig. 5, by the emission of several 
dislocations (seen as © Fig. 3(i)). The stress is relieved and a smoothing of the crack free 
surface, initially rough because of the void coalescent mechanism of growth (© in Figs. 3(d-i)), 
begins to take place resulting in a slight decrease of the surface area. Because the crack length is 
being measured through the crack free surface area, as described previously, this relaxation of 
the surface gives an effective “shortening” of the crack, producing a region of apparent negative 
growth after the crack length has reached its peak at t = 60 ps in Fig. 5. The distance between the 
two crack tips does not change in the process, indicating that the apparent “shortening” is not a 
result of crack closure, but only due to surface relaxation effects. 

Qualitatively, very similar behavior is observed for the prestress of 4.0 GPa (see the 
corresponding curve in Fig. 5). The average crack speed has increased to 701 m/s and the crack 
has stopped its growth at almost the same length of 35 mn (Fig. 4(b)), but earlier in time 
compared to the case of 3.75 GPa prestress. The region of the follow up negative growth for a = 
4.0 GPa is considerably deeper than for a = 3.75 GPa, indicating more intensive relaxation - a 
result of more intensive void formation and plastic processes accompanying the faster crack 
growth. 

When the prestress is increased up to 4.25 GPa another substantial change in the crack 
behavior is observed (see the corresponding curve in Fig. 5). Unlike simulations at lower values 
of prestress, the crack growth at o = 4.25 GPa occurs in three discrete increments. The first 
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increment occurs from 0 < t < 25 ps at an average speed of 752 m/s, while the second increment 
occurs from 75 < t < 100 ps at an average speed of 687 m/s, and the third at 125 <t< 155 ps at 
an average speed of 651 m/s. These repeating cycles of resumed growth can also be observed on 
the velocity vs. time curve shown in Fig. 6(a). 

A detailed analysis [34] revealed that the observed oscillatory behavior of the crack growth is 
a result of a dynamic instability [35,36] developed at the crack tip each time the instantaneous 
propagating velocity of the tip reaches a critical value v c of approximately 35% of the Rayleigh 
wave speed cr = 3180 m/s for this model. This instability was shown [34] to create a series of 
periodic dislocation bursts that appeared in the simulation shortly after the crack tip velocity 
reached v c . These dislocation bursts (some of which are still seen as © in Fig. 4(c)) reduce the 
crack propagation velocity. It is possible that the same mechanism of dislocation bursts is 
responsible for stopping the crack growth in the cases of a = 3.75 GPa and a = 4.0 GPa, where 
the tip velocity reached a peak of approximately 0.3 cr and 0.4 cr, respectively (Fig. 6(b,c)). The 
values of 0.3c« and 0.4c« may have been sufficient to induce dynamic instability and trigger the 
observed subsequent dislocation emission (© in Fig. 3(i) and Fig. 4(b), respectively), thus 
blunting the crack tip in the +x direction for both cases. In the case of 4.25 GPa prestress, the 
accumulated strain energy in the system was sufficient to overcome this dynamic blunting and 
restore the crack growth on two successive occasions, producing the velocity oscillations 
observed in Fig. 6(a). 

Another important factor in the process of fast intergranular crack growth is the role of the 
GB. When dynamic instability occurs, the GB interface prevents crack branching - a typical 
phenomenon in brittle fracture [35,36], making it possible for the crack to restore its growth 
when the propagation speed decreased below v c . Consequently, the crack can overcome the 
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dynamics instability and continue its growth if the system has a sufficient amount of strain 
energy. The case which is observed for the a = 4.25 GPa prestress (as seen in Fig. 4(c) and Fig. 
5). 

4. Energy release mechanisms during crack growth 

The process of crack growth, as described in Section 2, occurs in a prestrained system under 
fixed displacement boundary conditions. During crack growth, there is no external work done on 
the system and the system exchanges energy only in the form of heat transfer through the applied 
Nose-Hoover thermostat [20]. The crack grows by expending elastic energy that was stored 
when the system was initially prestrained. The energy balance can be expressed as 

^ E e, = E s ~ A E GB + E p ,+AQ, (7) 

where: A E el is the decrease of the stored elastic energy as the crack grows; E s = 2/y s , is the 
energy of the newly formed free surface of the crack; A E GB = /y G B, is the excess (i.e., above the 
perfect crystal) energy of the destroyed GB interface; E , is the energy related to all plastic 
processes taking place in the plastic zone around the crack tips; and A Q is the heat released in 
the process. 

By summing the energies from every atom in each of the structural states discussed in 
Section 2, it is possible to determine each structural state's contribution to the potential energy of 
the system (e.g., free surface energy, GB energy, elastic energy and plastic energy). The total 
free surface energy, E s in Eq. (7), is determined as the sum of the potential energy of all surface 
atoms (as defined in Section 2). The GB energy E CB is approximated as the sum of the potential 
energy of all disordered atoms. The computed change of the grain boundary energy, -A E GB in 
Eq. (7), is negative since the overall number of GB atoms decreases due to the transformation of 
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GB atoms into surface or crystalline atoms as the crack grows. Note that this calculation of grain 
boundary energy is approximate since a small number of disordered atoms are also present as 
either dislocation cores or vacancies. 

The elastic energy, E e i in Eq. (7), is determined as the excess energy of the atoms in a fee 
state, as identified by the CNA technique [29]. Here, the atoms in the fee state form the 
elastically deformed crystal lattice, and their excess potential energy is the elastic energy of the 
lattice. Plastic deformation occurs mainly through the formation of twins and stacking faults as a 
result of dislocation nucleation. Since twins and stacking faults are comprised of hep atoms [32], 
the plastic energy, E p i in Eq. (7), is determined as the sum of the excess potential energy of all 
the hep atoms in the system. However, this energy is only a part of the total energy consumed by 
the plastic processes. Another part, for example, the energy required to drive a dislocation, is 
transferred into heat and contributes to the A Q term in Eq. (7). 

The heat generated during crack growth is taken out of the system by the MD thermostat and 
by the absorbing layers as discussed in Section 2. The Nose-Hoover MD thermostat [20] 
introduces external non-conservative forces to the atoms decreasing or increasing their kinetic 
energy to keep the overall temperature constant. The work done by these non- conservative forces 
and by the damping forces applied to the atoms in the absorbing layers results in heat release of 
the system and can be calculated during the simulation to get the tenn A Q in Eq. (7). 

Despite its convenience and relatively easy implementation, the energy identification 
procedure outlined in the previous paragraphs, has some deficiencies. First; the precision of the 
identification depends ultimately on the precision of the CNA and the surface atoms 
identification techniques. The precision of both of these techniques depends on the ability to 
distinguish between the first and second nearest neighbors of a given atom. This distinction is 
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generally based on the distance between the neighboring atoms and is most appropriate when the 
temperature is relatively low (at room temperature and below when the thermal atomic vibrations 
are low), and when the elastic strain is relatively small (on the order of a few percent, so that the 
crystal lattice is not overly deformed). Second; the actual energetic contributions of each of the 
structural states is not as clearly delineated as this method suggests. For example, part of the 
potential energy of a non-fee atom should also be included in the elastic energy, as the non-fee 
structures also have elastic energy. Nevertheless, the energy to create such a structural defect in a 
perfect lattice (the excess energy) is usually much higher than the elastic energy of the 
formation*, thus latter can be neglected. Considering the issues outlined above, the suggested 
procedure in sorting out the different energy contributions in the process of crack growth is a 
good trade-off between simplicity, effectiveness and precision. 

Under the prestrained conditions of the model (Fig. 1), the stored elastic energy is the only 
source of energy for the crack to grow. The decrease of the elastic energy for a growing crack for 
each of the four applied prestress conditions is shown in Fig. 7. For comparison, FEM 
calculations for quasistatic crack growth in a pure elastic (no plasticity) material with the same 
elastic coefficients as those derived from the EAM interatomic potential [19] are given. In all 
cases, the MD energy curves fall below the FEM lines, indicating a faster release of the stored 
elastic energy than is expected for the quasistatic crack growth in the absence of plasticity. For 
the cases of 3.75 and 4.0 GPa prestress, the majority of energy release occurs at the end of the 
crack growth when the crack is blunted by a dislocation burst (see © in Fig. 3(i) and Fig. 4(b)). 
The process is accompanied by relaxation and smoothing of the initially rough crack surfaces, 


For example, the formation energy of a vacancy for the interatomic potential used in this study [19] is 0.64-^0.68 
eV/atom, while the strain energy for 1% strain is only 0.001 eV/atom. 
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when the crack becomes stationary, producing an effective shortening of the crack length as seen 
in Fig. 7 and Fig. 5. 

For the case of 4.25 GPa, there are three noticeable localized drops in the energy curve at 28 
mn, 45 mn and 60 nm. These instantaneous decreases in the stored elastic energy coincide with 
the time of the three successive peaks in the crack speed, shown in Fig. 6(a), when the critical 
velocity is reached. At this velocity, the fast propagating crack tip becomes dynamically unstable 
and produces dislocation bursts [34]. The observed (see Fig. 7) drops in the stored elastic energy 
for a = 4.25 GPa case are a direct result of the dynamically initiated dislocation bursts. 

As the crack grows, the elastic energy plus the energy of the destroyed GB is transferred into 
surface energy, plastic energy and heat (Eq. (7)). A plot of each of the energy contributions in 
Eq. (7), coming from different processes during crack growth are shown in Fig. 8. For the three 
given prestresses of a = 3.75, 4.0 and 4.25 GPa (in the case of a = 3.5 GPA the growth of the 
crack was insufficient to reliably determine the energy), the largest part of the energy was 
released in the form of heat. The next most substantial amount of energy was expended creating 
the free surface. A small amount of energy was recovered by destroying the GB. The potential 
energy of the hep atoms (given in Fig. 8) can be related to the stored plastic energy E pl in Eq. (7) 

as it gives the energy of formation for the twin boundaries and the stacking faults produced in the 
process of crack growth (seen in Figs 3 and 4). Fig. 8 shows that the energy associated with the 
hep atoms is small compared to the other forms of energy release. The greatest component of the 
energy released in the plastic processes does not go into the creation of these defects, but into 
viscous heat dissipation as defects are swept away from the crack tip. Another contribution to the 
heat release, shown in Fig. 8, comes from the phonons produced in the process of interface 
debonding. Since the heat release associated with phonons was not specifically determined, the 
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curves for the heat release in Fig. 8 reflect the total amount of heat of all dissipative processes in 
the system. 

The slope of the energy release curves, shown in Fig. 8, gives the energy release rate. Fig. 9 
compares the energy release rates of the heat release, free surface formation, and GB destruction 
for the four average initial crack tip velocities (estimated in Fig. 5) obtained at the four 
prestresses. As expected, the energy release rates for surface creation and GB destruction do not 
depend on the crack propagation speed. These rates (see Fig. 8) are close to twice the surface 
energy per unit area y s , and to the GB energy per unit area y GB , respectively (y s « 1 J/nr [19] and 
y gb ~ 0.65 J/nr [34]). In fact, y s extracted from Fig. 9 is slightly higher than 1.0 J/m 2 , which can 
be attributed to the non-relaxed state of the crack free surface. The GB energy released per unit 
area y GB also deviates slightly from the static value of 0.65 J/m 2 , and tends to decrease at high 
crack propagation speeds (above 700 m/s). As y GB should not depend on the speed of the crack 
propagation, the only conceivable explanation appears to be that the decrease of the disordered 
atoms, belonging to the GB, is partly compensated by the appearance of more disordered atoms 
via the formation of dislocations at high crack tip velocity (recall that the GB energy is defined 
as the energy of the disordered atoms, assuming that most of them belong to the GB interface). 

Unlike the energy release rate spent to create the crack free surface, the rate of heat release 
increases rapidly with increasing crack propagation speed (Fig. 9). While at the lowest simulated 
speed of 25 m/s, the heat release rate is lower than the surface energy release rate, at high speed 
it becomes the dominant contribution to energy release rate. Initially, Gheat increases linearly with 
the crack propagation speed, suggesting viscous crack propagation. At high velocities Gheat 
increases faster than linear, which is consistent with both the theoretical basis for the MD 


19 



Invited paper for the Special Issue of Journal of Materials Science 
Special Issue Title: “Nanostructured Materials and Mechanical Behavior” 

simulation model [37] and experimental observations of a crack moving in a brittle amorphous 
material [38], 

5. Conclusions 

Using a molecular-dynamics model for crack propagation under steady-state conditions, the 
presented paper studies the dynamics and energetics of intergranular crack growth along a flat 
grain boundary in aluminum. The crack is initiated by screening the interaction forces between 
atoms in a short region on both sides of the grain-boundary interface in a hydrostatically 
prestressed bicrystalline aluminum plate. The dynamics of the crack shows a strong non-linear 
dependence of the speed of the crack propagation on the prestress. When the crack propagation 
speed increases to 0.35% of the Rayleigh wave speed, a dynamic instability occurs at the crack 
tip producing periodic dislocation bursts until the crack is blunted. The GB interface prevents 
crack branching. 

Using the ability of the molecular-dynamics simulation to identify atoms in different 
structural configurations, it was possible to identify the energy contribution of different 
processes taking place during crack growth. The energy contributions were determined as: (i) 
elastic energy - defined as the excess potential energy of the atoms in fee crystallographic state; 
(ii) plastically stored energy - the excess energy of stacking faults and twin boundaries; (iii) 
grain-boundary energy and (iv) surface energy. Additionally, the energy dissipated as heat in the 
system was determined by monitoring the heat exchange through the thermostat. This energetic 
analysis indicates that the majority of energy in a fast growing crack is dissipated as heat. This 
dissipation increases in a linear fashion at low speed and faster than linear at a speed approaching 
1/3 the Rayleigh wave speed. At relatively low speed (below 200 m/s or less than few percent of 
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the Rayleigh wave speed), the dissipated heat becomes lower than the energy needed to create 
the crack free surface. 
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Figures: 



Fig. 1 . Schematic representation of the simulation model, explained in the text. Size dimensions 
are in mn. 
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Fig. 2. Atomistic snapshot giving the crystallography and structure of the GB interface. Common 
neighbor analysis (CNA) [29,30] is used to identify atoms in different crystallographic states: fee 
(small dots), hep (triangles), and non-crystalline atoms (large dots). Atoms with more than 1/3 of 
their nearest neighbors missing are identified as surface atoms (squares), indicating existing 
vacancies in the GB. The length scale is in units of the lattice constant of Al, a 0 = 0.405 nm. 
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Fig. 3. A series of snapshots monitoring the process of crack growth in the MD system 
prestressed at 3.75 GPa tensile hydrostatic load. The snapshots are taken at various times from 
the crack initiation, t, given in ps. CNA is used to identify atoms in hep crystallographic state 
(forming the grey lines of twin boundaries and stacking faults), and non- crystalline atoms (in 
black indicating the GB and dislocation cores). Thus, a number of different formations are 
indicated in the figure as follows: © - GB interface; © - deformation twin; (D - partial or 
twinning dislocation; © - nanovoid at the crack tip; © - slip dislocation. The wavy Moire 
patterns in the crystal regions make visible the strain patterns and phonon waves in the crystal 
(see text). 
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Fig. 4. MD snapshots of cracks, which have 
propagated in the MD system, shown in Fig. 1, 
prestressed at three different initial hydrostatic 
stresses: a = 3.5 (a), 4.0 (b) and 4.25 GPa (c). As 
in Fig. 3, CNA is used to identify atoms in 
different formations ©-©, which correspond to the 
types indicated in Fig. 3. 


27 



Invited paper for the Special Issue of Journal of Materials Science 
Special Issue Title: “Nanostructured Materials and Mechanical Behavior” 



Fig. 5. Crack length vs. time for the four prestress values applied in the simulation. 
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Fig. 6. The speed of crack growth in units of the Rayleigh wave speed, cr = 3180 m/s, vs. time 
for the four prestress values applied in the simulation. 
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Fig. 7. The stored elastic energy in the simulated system vs. the length of the growing crack. 
Symbols represent MD simulation results. The solid lines represent FEM quasistatic elastic 
simulation results. 
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Fig. 8. The energy released through four different mechanisms accompanying the crack growth: 
heat dissipation, free surface formation, growth of stacking faults and twin boundaries, and the 
energy gained by destroying the GB interface. 
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Fig. 9. The energy release rate for three different mechanisms accompanying the crack growth: 
heat dissipation, free surface formation, and the rate of energy gained by destroying the GB 
interface. 
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